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Abstract 



The implementation and reliability of a quadratic diffusion Monte Carlo 
method for the study of ground-state properties of atoms are discussed. We 
show in the simple yet non-trivial calculation of the binding energy of the 
Li atom that the method presented is effectively second-order in the time 
step. The fulfilment of the expected quadratic behavior relies on some ba- 
sic requirements of the trial wave function used for importance sampling, in 
the context of the fixed-node approximation. Expectation values of radial 
operators are calculated by means of a pure estimation based on the forward 
walking methodology. It is shown that accurate results without extrapolation 
errors can be obtained with a pure algorithm, explicitely reported, that can 
be easily implemented in any previous diffusion Monte Carlo program. 



I. INTRODUCTION 



Quantum Monte Carlo (QMC) methods have been widely employed in quantum chem- 
istry calculations |l|-f|]. Nowadays, QMC methods have achieved accuracies comparable 
with CI calculations, mainly for small atoms and molecules. QMC is specially well suited to 
obtain the ground-state energy, providing an accurate treatment of electronic correlations in 
atomic and molecular systems. Apart from some technical and/or algorithmic differences, 
all the QMC methods pursue a common objective which is to solve the Schrodinger equation 
by relating the Green's function of the system to a random walk. In the present work, we will 
focus on the diffusion Monte Carlo (DMC) method, based on a short time approximation to 
the Green's function. The required antisymmetry of the atomic and molecular ground-state 
wave functions leads to the well-known sign problem. The present analysis is restricted to 
the standard fixed-node (FN) approximation According to the FN method, a prefixed 
nodal surface is defined in the configuration space of the system. That model nodal surface, 
which does not change along the calculation, introduces a boundary condition that forbids 
any crossing through a node. It has been proved that the resulting FN energy is an upper 
bound to the exact ground-state energy J| , the quality of the upper bound being related to 
the overlap between the model and the exact nodal surfaces. 

The DMC method is simpler than domain Green's function Monte Carlo |J but, contrary 
to this, DMC presents a time-step At dependence that requires additional analysis. In gen- 
eral, the use of a finite At to approximate the imaginary-time Green's function introduces a 
bias in the calculation. That problem is addressed by using several small time steps and then 
extrapolating the results for the energy to the At — > limit. The simplest implementation 
of DMC uses a short-time Green's function which is accurate to order (At) 2 (Eq. |5|). In 
this case, one obtains energies that behave linearly with sufficiently small At and then an 
extrapolation to zero time-step is unavoidable. In order to reduce that large time-step de- 
pendence, and hence accelerate the rate of convergence of the calculation, several algorithms 
have been proposed |6|-|T0|. Reynolds et al. [[3],(4| introduced an accept/reject step after the 
gaussian and drift movements in order to exactly sample ($o) 2 for all time steps if the exact 
wave function <3> were used for importance sampling. Since in a true calculation the trial 
wave function is close to the exact one, the time-step dependence is significantly decreased 
with respect to the linear case. On the other hand, the rejection method introduces the 
problematic issue of persistent configurations , giving rise to a somehow unpredictable 
behaviour of the final At dependence of the energy. 

A second alternative to reduce the time-step bias is to work out second-order algorithms 
by considering short-time Green's function accurate to order (At) 3 (Eq. §). Some time ago, 
Chin || developed second-order propagators and proved the quadratic efficiency of them 
in several model problems. That approach, here referred to as quadratic DMC (QDMC) 
method, has been profusely used in the microscopic study of quantum fluids, mainly 4 He 
12|1 and 3 He In all these applications, QDMC has proved to be a real second-order 



algorithm and therefore the time-step bias is much more under control than in linear DMC 
algorithms. It is worth noticing that liquid 3 He is a fermion system, and the presence of 
a fixed nodal surface (FN approximation) has not broken down the expected second-order 



accuracy [13]. This feature contradicts some preliminary conclusions M that questioned 



a (At) 2 behavior when using nodal trial functions. Recently, Forbert and Chin |l4j] have 
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improved the time-step accuracy by developing a fourth-order algorithm that has manifested 
a fourth-order power law in a exigent calculation like the binding energy of liquid 4 He. 

The use of importance sampling (ip) in the DMC solution of the Schrodinger equation 
makes walkers to be distributed according to the mixed probability distribution ip$>o, and 
not ($o) 2 - Therefore, mixed estimators, which are the natural output in DMC, are only 
unbiased when the operators in question coincide with the Hamiltonian H or commute with 
it. In a first and very common approximation the bias of mixed estimators of, for instance, 
radial operators are dealt with an extrapolation method ||. Within this procedure, a less 
biased estimation is obtained from a proper linear combination of variational and mixed 
estimators. The efficiency of that linear extrapolation is strongly related to the variational 
quality of ■?/>, and the results so obtained are still biased in a quantity difficult to assess 
pT5| , [T6[1 . It is possible, however, to eliminate the uncertainties present in the extrapolation 
method using really pure estimators. 

The objective of obtaining pure estimators has been accomplished by using different 
techniques mainly bilinear sampling [ 17] , path integral formalism fl8| , |l9ll and forward walking 
p0| , p!5|j2T| , [22|JT6| , |23[| . Among them, several implementations of forward walking have been the 
most explored algorithms to achieve unbiased estimations of radial operators. Following 
Liu et al. |24|| , the necessary ingredient to correct the mixed estimator, i.e., the quotient 
($o(R-)/ V'(R'))) can be obtained from the asymptotic offspring of the R walker. In some 
implementations of that forward walking procedure, tagging algorithms p0|-f22|1 have been 
devised in order to estimate the asymptotic number of descendants. The complexity of an 
efficient tagging algorithm, together with the large fluctuations observed in the asymptotic 
offspring, have impeded its widespread use among the DMC community. However, there is a 
simpler algorithm to sample pure estimators which does not require any tagging algorithmic 
structure and that can be very easily incorporated in any DMC code [111. This method 



has been checked in model problems and has proved its accuracy in many calculations of 
quantum liquids properties. 

In this work, we address the two above discussed aspects of the DMC method: the 
problem of the time-step dependence and the sampling of the exact distribution to extract 
pure estimators. We propose the use of a second-order DMC algorithm, that presents 
small time-step errors, and unbiased estimators to calculate exact values for non differential 
properties. As it is commented in the next section, the implementation of that algorithm can 
be done rather straightforwardly in any DMC code. The reliability of the method has been 
analyzed by studying the ground state of atomic lithium. This simple but not trivial system 
has been chosen since, on one hand, presents all the ingredients to test the performance of 
the method, and on the other, there are very accurate calculations with which we 

can compare our results. 

The rest of the paper is organized as follows. In Sec. ||, the QDMC algorithm and the 
pure estimation methodology is presented. Sec |TTT| is devoted to present and discuss the 



results. Finally, the main conclusions of the present work are presented in Sec. IV 



II. METHOD 
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A. Quadratic Diffusion Monte Carlo 



The DMC method solves the Schrodinger equation in imaginary time for the function 
/(R,t) = V(R)$(R,t), 

- = ~\ V 2 R /(R, t)+ 1 - V R (F(R)/(R, t)) + (E L (R) - E) /(R, t) , (1) 

$(R, t) being the wave function of the system and ^(R) a trial function used for importance 
sampling. Equation ([I]) is written in atomic units, which are the ones used throughout this 
work. In Eq. (|), E L = ^(R) _1 ^(R) is the local energy and F(R) = 2^(R) _1 V R ^(R); 
R stands for a 3iV-coordinate vector and E is an arbitrary energy shift. 

The r.h.s. of Eq. ([!]) may be written as the action of three operators Oi acting on the 
wave function /(R, t), 

- = (d + 2 + O s ) / (R, t) = O /(R, t) (2) 

The three terms Oj are interpreted by similarity with classical differential equations. The 
first one, 0±, corresponds to a free diffusion with a diffusion coefficient D — 1/2; O2 acts 
as a driving force due to an external potential, and finally O3 looks like a birth/death 
term. In Monte Carlo, the Schrodinger equation (0) is best suited when it is written in a 
integral form by introducing the Green's function G(R',R, At), which gives the transition 
probability from an initial state R to a final one R' in a time At, 

/(R', t + At) = J G(R, R, At) /(R, t) dR . (3) 

More explicitly, the Green's function is given in terms of the operator O by 

G(R',R,At) = (R'| exp(-OAt) |R). (4) 

Most applications of the DMC method work with the simplest version of the short-time 
approximation, 

exp (-OAt) = exp (-0 3 At) exp (-0 2 At) exp {-O x At) + O ((At) 2 ) . (5) 

This expansion generates a linear time-step dependence that is reproduced in any calculation 
when the dependence of the energy on At is analyzed. A significantly better behavior is 
obtained by expanding the exponential of the operator O to higher orders in At. A good 
compromise between algorithmic complexity and efficiency is obtained by using a second- 
order expansion (QDMC). In this case, the Green's function G(R',R, At) is approximated 
by 

exp (-OAt) = exp (-0 3 -f \ ex P (~° 2 t) 6XP ( "° lAt) (6) 



x exp (-02^) exp {-Os^j + O ((At) 
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This decomposition, which is the one we have used in our calculations, is not unique as 
pointed out by Rothstein et al. [|7| and Chin || . Introducing the above expansion (6) in Eq. 
(H) the Schrodinger equation, written in integral form, becomes 



f(R',t + At) = / 



(/ :! j R', Ri, — \ G 2 ^Ri, R2, — ) ('\ (R-2- R > At) 



(7) 



x G 2 ( R3, R/i, G 3 (~R4, R 



2 J " V 2 7 



/(R,t) dRi...dR 4 dR . 



In Eq. (7), the total Green's function G is split into the product of individual Green's 
functions Gj, each one associated to the single operator 0{. G\ is the Green's function 
corresponding to the free diffusion term (Oi), and thus it is the well-known solution for a 
noninteracting system, 



:».V 



Gi(R',R,t) = (4nDt)-— exp 



(R/-R) 
4Dt 



(8) 



In the MC simulation, the evolution given by G\ corresponds to an isotropic gaussian move- 
ment of size proportional to \f~Dt. The Green's function G 2 describes the movement due to 
the drift force appearing in 2 ; its form is given by 

R(0) = R 

G 2 (R / ,R,t) =5(R'-R(t)), where <{ . (9) 



It 

Under the action of G 2 , the walkers evolve in a deterministic way according to the drift force 
F(R(t)). In order to preserve the second-order accuracy in the time step, the differential 
equation (H) must also be solved with a second-order integration method. Finally, the third 
individual Green's function G3 has an exponential form with an argument that depends on 
the difference between the local energy of a given walker and the prefixed value E, 

G 3 (R', R, t) = exp [-(E L (R) - E) t] 8(R' - R). (10) 

This third term, which is called the branching factor, assigns a weight to each walker ac- 
cording to its local energy. Depending on the value of this weight the walker is replicated or 
eliminated in the population list. A pseudo code which describes the evolution of a walker 
for a time step At is presented in Ref. [|29[ . It is worth noticing that a second-order algo- 
rithm, like the one presented here, does not introduce any new terms with respect to a linear 
approximation. Therefore, it is quite simple to move from a linear code to a second-order 
one from the programmer's point of view. 

As mentioned earlier, instead of using higher-order evolution operators some authors 
have been using a linear DMC method which incorporates an accept/reject step (Metropolis 
DMC -MDMC-) After a proposed gaussian plus drift movement for an individual 

electron, the new position is accepted with probability 

. ( \i.(R')\'G(R,R',At) \ 
p = nm U(R)|»G(R',R,At)' 1 ) ' ( ) 
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The detailed balance condition (|TT|) ensures theoretically a more accurate sampling of \ip\ 2 , 
but introduces in DMC the potential problem of persistent configurations. Results obtained 
with the MDMC method show a significant reduction of the time-step bias, with respect to 
DMC, but its efficiency depends on the particular system under study and the extrapolation 
law to At = can show erratic behaviors pO . 



B. Pure Estimators 

When the asymptotic limit (t — > oo) is reached, the sampling of an operator O is carried 
out according to the mixed distribution V^o, with $ the ground-state wave function. Thus, 
the natural output in DMC corresponds to the so called mixed estimators. The mixed 
estimator of an operator O(R) is, in general, biased by the trial wave function ip used 
for importance sampling. Only when O(R) is either the Hamiltonian, or commutes with 
it, the mixed estimation is the exact one. A simple method that has been widely used to 
approximately remove the bias present in the mixed estimations is the extrapolated estimator 



(0(R)) e = 2(0(R)) m -(0(R)) e , 
from the knowledge of the mixed estimator (0(R)) m and the variational one 

(0(R)>„ 





O(R) 


^(R)> 


MR) 


V(R)> 



(12) 



(13) 



The expectation values obtained through the extrapolation method (|T2| ) depend on the 
trial wave function ip used for importance sampling. Therefore, in spite of using good trial 
wave functions, the extrapolated estimator is always biased in a quantity difficult to assess 
a priori. In order to overcome that important restriction, one can go a step further and 
calculate pure (exact) expectation values, 



(0(R)>. 



($o(R) | O(R) | $ (R)) 
($o(R) | $ (R)> 



Having in mind that walkers evolve according to the mixed distribution 
estimator is more conveniently written as 



(0(R)) P = ( $„(R) 



O(R) 



$o(R) 
V(R) 




$o(R) 



go(R) 
^(R) 



V(R) 



(14) 
the pure 

(15) 



Some time ago, Liu et al. [24] proved that $o(R)/V ; (R) can be obtained from the asymp- 
totic offspring of the R walker. Assigning to each walker Rj a weight W(Rj) proportional 
to its number of future descendants 



Eq. flT5p becomes 



W(R) = n(R,t-> oo) , 
EiO(Ri) W(R< 



(0(R)) P 



(16) 



(17) 
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where the summatory Ya runs over all walkers and all times in the asymptotic regime. The 
difficulty of the method, known as forward walking, lies on the estimation of the weight 
VT(R) (|16|). The weight of a walker existing at time t is not known until a future time 
t' > t + T, T being a time interval long enough so that Eq. fll6"f) could be replaced by 
W(R(t)) =n(R(t')). 

The evaluation of Eq. (|T7|) has traditionally required the implementation of a tagging 
algorithm pp|-fm. The purpose of that algorithm is to know, at any time during the sim- 



ulation, which walker of any precedent configuration originated a present walker. In this 
way, one can determine the number of descendants of each R i; and accumulate its contri- 
bution to Eq. fllTD from the distance. An easier method consists in working only with the 
present values of O(Rj), in such a way that a weight proportional to its future progeny is 
automatically introduced E6|. The schedule of this second algorithm is the following. The 
set of walkers at a given time {R«}, and the values that the operator O takes on them {Oj}, 
evolve after a time step to {R^} and {O^}, respectively. In the same time step, the number 
of walkers N changes to N'. In order to sample the pure estimator of O, we introduce an 
auxiliary variable {Pi}, associated to each walker, and that evolves with it (its final average 
value will provide the pure expectation value of O). The evolution law for {Pi} is given by 

{p^ - {p;} = {p-} + m , (is) 

where {Pf} is the old set {Pi} transported to the new one, in the sense that each element Pi 
is replicated as many times as the Rj walker, without any other changes. {Pi} is initialized 
to zero when the run starts. 

In order to better fulfill the asymptotic condition ([[(]) the evolution of the {Pi} values 
is further carried on by means of 

{Pi} - {Pi} = {Pf} • (19) 

Since a calculation is usually divided in blocks, one can collect/transport data during a block 
via Eq. (|T8|), and allow for a only reweighting by means of Eq. ([19]) in the following one. 
Thus, after a first initialization block, each new block gives a value for the pure expectation 
value of O 

N f 

(0(R)>„ = £{P<}/(MxJV» , (20) 
i-i 

M being the length of the block and Nr the number of walkers at the end of each reweighting 



block. A pseudo code for the pure estimation is available in Ref . |29| . It is worth mentioning 
that the algorithm can be easily generalized to DMC programs that work with branching 
weights. 

The pure estimation depends on the size M of the block. M has to be large enough to 



match the asymptotic regime required by the forward walking property (16). With an easy 
modification of the code reported in Ref. , one can introduce a pure estimator for a set of 
increasing values of M in a single calculation. Looking at the results obtained, as a function 
of the block size, one observes a characteristic value M c from which on the pure estimator 
is the same (considering the respective statistical errors). 
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III. RESULTS 



The efficiency and accuracy of both the QDMC method and the pure estimator described 
in the previous section are analyzed in a calculation of ground-state properties of atomic 
Li. This is a well-known system, with nearly exact results, and serves well to our purpose 
since it contains the ingredients that can hinder the pursued second-order behavior: a nodal 
surface and cusp conditions. The present calculation works in the fixed-node approximation 
and therefore the energy obtained is an upper-bound to the ground-state energy of the atom. 
The trial wave function used for the importance sampling is the usual Jastrow-Slater form 

V^j^, (21) 

with 0^ (0^) the Slater determinant of spin up (spin down) electrons. The two-body Jastrow 
factor is taken as 

*-g-(r^;)- (22) 

with coefficients CV,- = 1/2 (CV,- = 1/4) for spin- unlike (spin-like) electrons to satisfy the 
electron-electron cusp conditions. 



The functional forms for the orbitals entering in are taken from Ref. [p3|| . 

^ls _ g-Ci s r e -wr/{l+vr) ^3) 

x 2s = (1 + C r) e~^ r e -W(i+w) . 

They introduce an additional exponential, with respect to hydrogenic-like orbitals, to better 
fit the electron-nuclear correlations. The fulfilment of the electron-nuclear cusp condition 
introduces two constraints that reduce the number of parameters to be variationally opti- 
mized to four: b, v, and ( 2s - We have taken for them the values reported by Langfelder 



et al. 23 



The VMC energy obtained with i[) is in agreement with the one reported in Ref. |23 



E = —7.4737(1). That supposes a 90% of the correlation energy, which is significantly high 
taking into account the relatively simple trial wave function used. That result coincides 



with the more sophisticated wave function vpf proposed by Schmidt and Moskowitz [3T 
and is slightly worse than the energy obtained using (E = —7.47669) from the same 
authors. In a DMC calculation, like the one intended here, the proposal of Langfelder et 



al. [^] is more appropriate since it satisfies both the electron-electron and electron-nuclear 
cusp conditions. These requirements are essential in order to eliminate divergences of the 
drift force that break down the expected time-step accuracy. In fact, we have verified that 
if ^7 is used instead of ip a second-order behavior is not achieved. 

One of the main concerns of the present work is the time-step dependence of the different 
DMC algorithms discussed in Sec. II. In Table I, results for the energy as a function of At 
are reported using linear DMC, linear DMC with rejection (MDMC) and second-order DMC 
(QDMC). The exact non relativistic ground state energy is -7.478 060 323 10(31), calculated 
variationally by using a multiple basis set in Hylleraas coordinates |]26|| . The last row in 
Table I contains the extrapolation to At = based on a polynomial fit to the data. The 
QDMC results follow the expected second-order behavior 



S 



E(At) = E + E 2 (At) 2 (24) 

with good accuracy (x 2 /V = 1.2). The extrapolated energy E = —7.47805(2), which 
reproduces the exact value, is statistically indistinguishable from the energy obtained with 
At = 0.003. In fact, one may move in a range At = 0.0025 — 0.0035 without significant 
changes in the energy. This is a very useful feature which derives from the second-order 
Green's function used in the QDMC method, and that has been also observed in quantum 
liquids calculations [|T2],|T3]]. 

The extrapolations corresponding to the DMC and MDMC calculations have been carried 
out with a second-degree polynomial that, differently from the QDMC case (24]), includes a 



linear term. The variance of the DMC result is larger than the QDMC one and the energy 
obtained is slightly above the QDMC and exact energies. The different dependence on At 
of the DMC and QDMC algorithms is illustrated in Fig. la: DMC energies are essentially 
linear while QDMC ones are quadratic. 

The introduction of an accept/reject step in the imaginary-time evolution modifies com- 
pletely the time-step dependence of the DMC method. As is clear from Table I and Fig. 
lb the MDMC method allows for a significant reduction of the time-step bias, mainly for 
large At. For values At > 0.01, the MDMC energies are closer to the exact result than 
the QDMC ones. However, when At is reduced the behavior of the energy shows a quite 
erratic behavior (see Fig. lc) that makes difficult the extrapolation to zero time-step. The 
extrapolated value reported in Table I comes from a second-order polynomial fit. The value 
obtained, -7.47810(14) is in agreement with the QDMC result but its variance is an order 
of magnitude larger. 

The pure estimator described in Sec. II, in conjunction with the QDMC method, has 
been used to calculate expectation values of several radial operators. In particular, we have 
considered the first radial moments of both the one-body density and the two-body intracule 



and extracule densities 321, defined as 



(r«> 
(r5> 



N 

i=l 

N 

(*| E r£|*) 

i<j=l 

(*|*) 



N 

E_ R n \*) 

(R n ) = . (25) 



In Eqs. (|25|) , is the position vector of the z-electron, = — r,, and R = (r, + Tj)/2 
is the center of mass vector of electrons % and j. Those observables play an important role 
in the description of the structure and dynamics of atoms and molecules p2|,p3| . They have 



been calculated within different theoretical frameworks [p5| - |28|j3^ -p4] to determine electronic 
properties, to interpret the Hund rules, to elucidate electron correlation effects, or to study 
elastic and inelastic form factors. 

In Fig §, we show the dependence of the pure expectation values of (r), (r 12 ), and (R) 
upon the number of time steps M of a block. In these calculations, more CPU time has been 
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invested for the largest values of the block length M in order to reduce the statistical noise. 
Nevertheless, the pure estimators reach very soon the asymptotic regime, for M values as 
low as M ~ 500. The statistical error of the pure estimators, evaluated in the asymptotic 
regime, is typically less than twice the statistical error of the mixed estimators for a fixed 
CPU time. The asymptotic results, i.e., the real pure estimators, reproduce accurately the 
exact values. As a matter of comparison, the mixed estimations and the variational ones are 
also shown. With the mixed and variational estimators, one easily obtains the extrapolated 
values (12) that are also plotted in Fig. 2. In all cases, the extrapolated estimator does not 
match the exact values, manifesting a bias that depends on the particular operator. The 
above features are also applicable to the second-degree radial operators (r 2 ), (r 2 2 ), and (R 2 ) 
which appear reported in Fig. 3. 

Our best pure estimates for the (r n ), (r™ 2 ), and (R n ) moments are contained in Tables II, 

III, and IV, respectively. Within the statistical errors, all the pure results are in agreement 
with the values considered as exact. With nearly the same variance than the extrapolated 
values, pure estimators eliminate the bias which comes from the trial wave function. On 
the other hand, extrapolated values retain some bias from ip which is certainly difficult to 
assert a priori. It is worthwhile determining the accuracy of pure estimators using different 
propagators like DMC or MDMC. For that purpose, we have also implemented the pure 
estimator algorithm to the MDMC program. The results obtained are reported in Tables II- 

IV. The pure MDMC results improve the extrapolated values but do not share the excellent 
quality obtained with the QDMC algorithm. Somehow, the reason for that is the use of 
effective time-steps in the branching factor, the accuracy of which is a crucial point to 
ensure the success of the method. 



IV. CONCLUSIONS 

In the present work, the time-step dependence of the diffusion Monte Carlo method 
and the implementation of pure estimators, both in atomic systems, have been addressed. 
The study has been carried out for Li since it already contains the necessary ingredients 
to test the method and, additionally, a lot of nearly exact data are available. Using a 
short-time Green's function accurate to O (At) 3 , we have proved that it is possible in atomic 
systems to have a quadratic behavior in the simulated energy as was reported in quantum- 
liquid calculations. To the best of our knowledge it is the first stringent test of QDMC on 
atomic systems with nodal surface. The accuracy obtained had been occasionally questioned 
due to the presence of the divergences that appear in the drift force near the nodes. The 
present results show that the presence of nodes is not a real handicap, a conclusion that 
coincides with the one achieved in a recent QDMC calculation of homogeneous (fermion) 



liquid He |13 |. The most crucial point to guarantee a success of the QDMC method is 



the correct fulfilment of the electron-electron and electron-nuclear cusps conditions. The 



energies obtained with the relatively simple model for ip, proposed by Langfelder et al. p3 
follow accurately the expected second-order behavior. Additional checks performed using 
better trial wave functions, from the variational view, but not satisfying the electron-nuclear 
cusp show a more erratic behavior of the function E(At). 

The second aim of the work has been the consideration of a rather simple pure estimator, 
based on the forward walking methodology. Pure results for different radial moments are in 
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agreement with exact values and are manifestly superior to the extrapolated estimations. 

As a concluding remark, we think that the combination of a QDMC method and the 
pure estimation here presented constitutes an appropriate tool for studying atoms and small 
molecules. Moreover, and very importantly from the programmer's viewpoint, the effort i) 
moving from a usual DMC to QDMC, and ii) introducing a pure estimator, is really small. 
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TABLES 



TABLE I. Energy of the Li atom versus the time step and using different algorithms. The 
figures enclosed in parenthesis are the statistical errors. The exact energy is -7.4780632310(31) 



At PMC MDMC QDMC 

0.030 -7.45279(4) -7.47860(4) -7.48244(3) 

0.0200 -7.46074(5) -7.47844(3) -7.48007(3) 

0.0100 -7.46893(5) -7.47831(4) -7.47855(4) 

0.0075 -7.47127(4) -7.47810(3) -7.47831(3) 

0.0050 -7.47333(3) -7.47811(3) -7.47817(2) 

0.0030 -7.47520(4) -7.47827(4) -7.47807(4) 

extrapolation -7.47784(11) -7.47810(14) -7.47805(2) 
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TABLE II. Estimations of (r n ) moments of the Li atom using the MDMC and QDMC algo- 



rithms. 

(r) (r 2 ) (r 3 ) 

Variational* 4.9553(7) 18.088(7) 91.12(6) 

Extrapolated (MDMC) 4.992(3) 18.39(3) 92.8(3) 

Extrapolated (QDMC) 5.001(3) 18.46(3) 93.3(3) 

Pure (MDMC) 4.985(2) 18.32(2) 92.4(2) 

Pure (QDMC) b 4.991(2) 18.36(2) 92.7(2) 

exact 4.989523 c 18.354615 c 92.60364 d 



a VMC results of Ref. ||: (r) = 
b GSD-DMC results of Ref. ||: 
c Hylleraas expansion of Ref. [26] 
d Hylleraas expansion of Ref. |2 



4.9536(6), (r 2 ) = 18.076(6) 
(r) = 4.992(7), (r 2 ) = 18.36(7) 
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TABLE III. Estimations of (rf 2 ) moments of the Li atom using the MDMC and QDMC algo- 
rithms. 



(rh) 



(rf 2 ) 



Variational 
Extrapolated (MDMC) 
Extrapolated (QDMC) 
Pure (MDMC) 
Pure (QDMC) 
exact 



8.593(2) 
8.673(7) 
8.686(7) 
8.662(3) 
8.671(3) 
8.668397 a 



36.32(1) 
36.89(6) 
37.02(6) 
36.79(3) 
36.87(3) 
36.847838 a 



189.2(1) 
192.4(5) 
193.5(5) 
191.6(3) 
192.4(3) 
192.10037 b 



a Hylleraas expansion of Ref. [26 
b Hylleraas expansion of Ref. 
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TABLE IV. Estimations of (R n ) moments of the Li atom using the MDMC and QDMC algo- 
rithms. 



<J2> (&) (R 3 ) 

Variational 4.2615(7) 9.009(3) 23.48(2) 

Extrapolated (MDMC) 4.301(3) 9.15(2) 23.9(4) 

Extrapolated (QDMC) 4.309(4) 9.19(2) 24.05(7) 

Pure (MDMC) 4.309(3) 9.19(2) 23.76(7) 

Pure (QDMC) 4.300(2) 9.150(8) 23.89(4) 

exact 4.2996(6) a 9.145(3) a 23.86(1) 

l Ref. H 
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FIGURES 



FIG. 1. Total energy of the Li atom versus the time step At computed by using the DMC, 
MDMC and QDMC algorithms. The lines on top of the data correspond to polynomial fits: 
E(At) = E + Ex At + E 2 (At) 2 for DMC and MDMC, and E(At) = E + E 2 (At) 2 for QDMC. In 
all cases, the error bars are smaller than the size of the symbols. The dotted line stands for the 
exact energy. In a) the time-step dependence of the DMC and QDMC energies are compared. In 
b) the comparison is between QDMC and MDMC. Finally, c) illustrates the different behavior of 
MDMC and QDMC when At — >■ 0; notice the dispersion of the MDMC data compared with the 
regular behavior of the QDMC energies. 
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FIG. 2. Pure expectation values of first-degree radial moments of the Li atom as a function 
of the block length M. The solid, dashed, and dotted lines stand for the exact, variational, and 
extrapolated results, respectively. 
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FIG. 3. Pure expectation values of second-degree radial moments of the Li atom as a function 
of the block length M. The solid, dashed, and dotted lines stand for the exact, variational, and 
extrapolated results, respectively. 
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